"""
核反应堆子通道分析程序的物理计算函数
"""

import numpy as np
from utils.helper_functions import (
    get_fluid_properties,
    get_saturation_temperature,
    get_latent_heat,
    calculate_hydraulic_diameter,
    calculate_reynolds_number,
    calculate_friction_factor,
    calculate_nusselt_number
)
from constants.constants import (
    GRAVITY,
    SUBCOOLED_BOILING_COEFF,
    TURBULENT_MIXING_COEFF,
    INTERFACE_DRAG_COEFF
)

def calculate_phase_change_rate(pressure, T_l, T_v, T_w, alpha_l, alpha_v, 
                              channel_type, heated_perimeter, flow_area):
    """
    计算相变率
    
    参数:
    pressure: float, 压力 (Pa)
    T_l: float, 液相温度 (K)
    T_v: float, 气相温度 (K)
    T_w: float, 壁面温度 (K)
    alpha_l: float, 液相分数
    alpha_v: float, 气相分数
    channel_type: str, 子通道类型
    heated_perimeter: float, 加热周长 (m)
    flow_area: float, 流通面积 (m²)
    
    返回:
    float: 相变率 (kg/m³·s)
    """
    try:
        # 获取饱和温度
        T_sat = get_saturation_temperature(pressure)
        
        # 获取物性值
        props_l = get_fluid_properties(pressure, T_l, 'liquid')
        props_v = get_fluid_properties(pressure, T_v, 'vapor')
        
        if props_l is None or props_v is None:
            return 0.0
            
        # 计算潜热
        h_fg = get_latent_heat(pressure)
        
        # 根据子通道类型调整相变系数
        if channel_type == 'corner':
            boiling_factor = 0.8
        elif channel_type == 'edge':
            boiling_factor = 0.9
        else:  # center
            boiling_factor = 1.0
            
        # 计算相变率
        if T_w > T_sat:  # 沸腾
            # 计算沸腾换热系数
            if T_l < T_sat:  # 过冷沸腾
                h = calculate_subcooled_boiling_htc(T_w, T_sat, T_l, pressure)
            else:  # 饱和核态沸腾
                h = calculate_nucleate_boiling_htc(T_w, T_sat, pressure)
                
            q_boiling = h * (T_w - T_sat) * heated_perimeter
            gamma = boiling_factor * q_boiling / (h_fg * flow_area)
            
        elif T_v < T_sat:  # 凝结
            h = calculate_condensation_htc(T_v, T_sat, pressure)
            q_condensation = h * (T_sat - T_v) * heated_perimeter
            gamma = -boiling_factor * q_condensation / (h_fg * flow_area)
            
        else:
            gamma = 0.0
            
        return gamma
        
    except Exception as e:
        print(f"计算相变率时出错: {str(e)}")
        return 0.0

def calculate_subcooled_boiling_htc(T_w, T_sat, T_l, pressure):
    """
    计算过冷核态沸腾换热系数
    
    参数:
    T_w: float, 壁面温度 (K)
    T_sat: float, 饱和温度 (K)
    T_l: float, 液相温度 (K)
    pressure: float, 压力 (Pa)
    
    返回:
    float: 换热系数 (W/m²·K)
    """
    # 使用Chen相关式计算过冷沸腾换热系数
    h = SUBCOOLED_BOILING_COEFF * np.exp(2*pressure/1260 * (T_w-T_sat)/(2722*(T_w-T_l)))
    return h

def calculate_nucleate_boiling_htc(T_w, T_sat, pressure):
    """
    计算饱和核态沸腾换热系数
    
    参数:
    T_w: float, 壁面温度 (K)
    T_sat: float, 饱和温度 (K)
    pressure: float, 压力 (Pa)
    
    返回:
    float: 换热系数 (W/m²·K)
    """
    # 使用简化的核态沸腾相关式
    h = np.exp(2*pressure/1260 * (T_w-T_sat)/2722)
    return h

def calculate_condensation_htc(T_v, T_sat, pressure):
    """
    计算凝结换热系数
    
    参数:
    T_v: float, 气相温度 (K)
    T_sat: float, 饱和温度 (K)
    pressure: float, 压力 (Pa)
    
    返回:
    float: 换热系数 (W/m��·K)
    """
    # 使用Nusselt凝结膜理论
    props_l = get_fluid_properties(pressure, T_sat, 'liquid')
    props_v = get_fluid_properties(pressure, T_sat, 'vapor')
    
    if props_l is None or props_v is None:
        return 0.0
        
    rho_l = props_l['density']
    k_l = props_l['thermal_conductivity']
    mu_l = props_l['viscosity']
    h_fg = get_latent_heat(pressure)
    
    h = 0.943 * (rho_l * (rho_l - props_v['density']) * GRAVITY * k_l**3 * h_fg / 
                 (mu_l * (T_sat - T_v)))**0.25
    return h

def calculate_single_phase_htc(velocity, hydraulic_diameter, pressure, temperature, phase):
    """
    计算单相对流换热系数
    
    参数:
    velocity: float, 流速 (m/s)
    hydraulic_diameter: float, 水力直径 (m)
    pressure: float, 压力 (Pa)
    temperature: float, 温度 (K)
    phase: str, 相态 ('liquid' 或 'vapor')
    
    返回:
    float: 换热系数 (W/m²·K)
    """
    # 获取流体物性
    props = get_fluid_properties(pressure, temperature, phase)
    if props is None:
        return 0.0
        
    # 计算无量纲参数
    Re = calculate_reynolds_number(velocity, hydraulic_diameter, 
                                 props['density'], props['viscosity'])
    Pr = props['viscosity'] * props['specific_heat'] / props['thermal_conductivity']
    
    # 计算努塞尔数
    Nu = calculate_nusselt_number(Re, Pr)
    
    # 计算换热系数
    h = Nu * props['thermal_conductivity'] / hydraulic_diameter
    return h

def calculate_wall_friction(velocity, hydraulic_diameter, pressure, temperature, phase):
    """
    计算壁面摩擦力
    
    参数:
    velocity: float, 流速 (m/s)
    hydraulic_diameter: float, 水力直径 (m)
    pressure: float, 压力 (Pa)
    temperature: float, 温度 (K)
    phase: str, 相态 ('liquid' 或 'vapor')
    
    返回:
    float: 壁面摩擦力 (N/m³)
    """
    # 获取流体物性
    props = get_fluid_properties(pressure, temperature, phase)
    if props is None:
        return 0.0
        
    # 计算雷诺数
    Re = calculate_reynolds_number(velocity, hydraulic_diameter, 
                                 props['density'], props['viscosity'])
    
    # 计算摩擦因子
    f = calculate_friction_factor(Re)
    
    # 计算壁面摩擦力
    tau_w = 0.5 * f * props['density'] * velocity * abs(velocity)
    
    return tau_w * 4 / hydraulic_diameter

def calculate_interfacial_force(alpha_l, rho_l, u_l, u_v, v_l, v_v, hydraulic_diameter):
    """
    计算界面力
    
    参数:
    alpha_l: float, 液相分数
    rho_l: float, 液相密度 (kg/m³)
    u_l, u_v: float, 液相和气相轴向速度 (m/s)
    v_l, v_v: float, 液相和气相横向速度 (m/s)
    hydraulic_diameter: float, 水力直径 (m)
    
    返回:
    float: 界面力 (N/m³)
    """
    # 计算相对速度
    u_r = u_v - u_l
    v_r = v_v - v_l
    V_r = np.sqrt(u_r**2 + v_r**2)
    
    # 计算界面面积密度
    d_b = 0.001  # 气泡直径 (m)
    a_i = 6 * (1-alpha_l) / d_b
    
    # 计算界面力
    F_i = 0.5 * INTERFACE_DRAG_COEFF * rho_l * V_r * V_r * a_i
    
    return F_i

def calculate_turbulent_mixing(alpha_i, alpha_j, rho, dy, channel_type):
    """
    计算湍流混合
    
    参数:
    alpha_i, alpha_j: float, 相邻子通道的相分数
    rho: float, 密度 (kg/m³)
    dy: float, 横向网格间距 (m)
    channel_type: str, 子通道类型
    
    返回:
    float: 湍流混合质量流量 (kg/m²·s)
    """
    # 根据子通道类型调整混合系数
    if channel_type == 'corner':
        mixing_factor = 0.7
    elif channel_type == 'edge':
        mixing_factor = 0.85
    else:  # center
        mixing_factor = 1.0
        
    # 计算湍流混合
    L_tur = dy  # 湍流混合长度
    alpha_eq = 0.5  # 平衡空泡分数
    
    w_prime = (TURBULENT_MIXING_COEFF * mixing_factor * dy * L_tur * rho * 
              (alpha_j - alpha_i - (alpha_j - alpha_eq)))
    
    return w_prime
